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Abstract 

We present an R/C implementation of optimal simultaneous diagonalization of 
several real symmetric matrices using Jacobi plane rotations, with compact triangular 
storage of symmetric matrices. 
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Note: This is a working paper which will be expanded/updated frequently. All suggestions 
for improvement are welcome. The directory deleeuwpdx.net/pubfolders/sjacobi has a pdf 
version, the bib hie, the R and C code, and the complete Rmd hie. 


1 Introduction 

Suppose Aj are m real symmetric matrices of order n. Our problem is to hnd a square 
orthonormal K, i.e. a matrix with KK' = K'K = /, such that the matrices Hj = K'AjK 
are as diagonal as possible in the least squares sense. Thus we want to minimize the sum of 
squares of all off-diagonal elements of the Hj, or, equivalently, maximize the sum of squares 
of all diagonal elements. 

It is a natural idea to use Jacabi plain rotations, the orthonormal version of cyclic coordinate 
descent, for this. The algorithm was proposed, possibly first, by De Leeuw and Pruzansky 
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(1978). It was subsequently programmed in FORTRAN for the Applied Statistics algorithms 
by Clarkson (1988). Numerical analysts arrived later on the scene, but contributed a number 
of essential refinements and extensions, such as an algorithm for Hermitian matrices and a 
proof of quadratic convergence (Bunse-Gerstner, Byers, and Mehrmann (1993), Cardoso and 
Souloumiac (1996)). There is a version of the algorithm in R (R Core Team (2018)) in De 
Leeuw and Ferrari (2008), together with a number of related algorithms using Jacobi plane 
rotations. De Leeuw (2017) has an R/C program (i.e. calling routine in in R, computational 
routines in C) for m — 1, the classical Jacobi algorithm, but using compact triangular storage 
of symmetric matrices. In this paper we provide an R/C compact storage version of the 
simultaneous diagonalization algorithm. 

2 Algorithm 

In the classical Jacobi method we cycle through all off-diagonal elements in a fixed order, 
minimizing the sum of squares over single-parameter plane rotations each time. After a cycle 
is completed we test for convergence. If there is still room for improvement we start a new 
cycle. 

The plane rotation K for any i 7 ^ j is defined by 


v 


kjj — x i 
kji Vi 


where x 2 + y 2 = 1. For i ^ i, j we have ku— 1, and for z/ 7 - i,j and i 7 ^ i, j we have k v e — 0. 
For H = K'AK we find 


(z (,.r 2 ciijxy T o^jjy , 

(ijj (x y ) J- (da cijj)xy, 

objjX 2 T 2ciijXy T (i 11 y 2 . 

which implies 

hu + h 2 jj + 2 h 2 j = a| + a 2 j + 2a 2 J . 

For £ 7 - i,j 


ha 


hij hji 

hjj = 


ha hgi (1 /{X dj^y, 
hje h^j aay T djix. 


which implies 

1 2 1 ; 2 2 1 2 

Ki + hj t — a u + a je . 
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Thus for m of such matrices A k we want to find the optimal plane rotation by minimizing 


f(0) = ( a ijk(x 2 - y 2 ) + (o*ifc - a jjk)xyf 


k =1 


( 1 ) 


over the circle x 2 + y 2 = 1. Now define 


u = x 2 — y 2 , 
v = 2xy. 


( 2 ) 

( 3 ) 


Note that w 2 + u 2 = ( x 2 + y 2 ) 2 and thus u 2 + u 2 = 1 if and only if x 2 + y 2 = 1. In other words 
the map 


X 


u 




_y_ 


V 


maps a point on the circle to another point on the circle. If u 2 + v 2 = 1 then the inverse map 
consists of two antipodal points on the circle, of which one is 


x = 


'l + u 
2 ’ 


y = sign(uj 


1 — ii 


( 4 ) 

( 5 ) 


By (1) both antipodal points give the same loss function value, so we always choose the one 
in (4) and (4). 

With our new variables the problem of finding the optimal rotation is 


f(u,v) = J2( ua ijk + vd ijk f 


k =1 


with dij k = 4(djj k — a jjk)- Define the 2x2 matrix S by 


5 = 


p q 

q r 


( 6 ) 


with 
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P 'y y ®ijk > 

fc=l 

m 

Q ^ ' Q’ijkdijki 
k= 1 
m 

r = J2 d lk- 

k =i 

Then 

f(u, v) = pu 2 + 2 quv + ru 2 , (7) 

which we must minimize over w 2 + u 2 = 1. The minimizing (w,u) is a normalized eigenvector 
corresponding with the smallest eigenvalue of S. 

The two eigenvalues of S are 


A m ax(S) = ^{(p + r) + ^(p~r) 2 + 4q 2 }, 

Amin (S) = ^{(> + r) - \J(p- r) 2 + 4q 2 }. 

There are several cases to consider, depending on the precise structure of S. 

• If q 7 ^ 0 the two eigenvalues of S are different, and A m i n (S') < min(p, r). One choice for 
the corresponding normalized eigenvector has elements 


_ q _ 

\Jq 2 + (Amin(S') - P ) 2 

Amin ('S') P 

\Jq 2 + (A m in(5*) -P) 2 


( 8 ) 

(9) 


For this choice both u and v are non-zero, with v < 0. The eigenvector is unique, except 
for sign changes, which do not influence the loss function value in (7). So we always 
choose v < 0. From (4) and (5) we see that we can choose the solution 


K 


y2(l + u) - 1 / 2(1 -u) 

O 2 * 1 -") li/ah+T 


which has f(u,v ) = A m i n (5')- 


( 10 ) 
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• If q = 0 and p > r then u = 0 and v = ±1. Thus f(u,v ) = A min (S) = r, and it does 
not matter which sign we choose. We set 



• If q = 0 and p < r then v = 0 and u = ±1. We choose u = 1 in (10) , which makes K 
equal to the identity, and which means we skip this rotation. 

• If q = 0 and p — r then f(u, v ) = p + r for any choice of u and v. Again we set K equal 
to the identity. 

3 Illustration 

Suppose we have three 2x2 symmetric matrices 

## C, 1] [, 2] 

## [1,] +1 -1 

## [ 2 ,] -1 +1 

## C, 1] [,2] 

## [ 1 ,] +2 +0 

## [ 2 ,] +0 +0 

## C, 1] [, 2] 

## [ 1 ,] +1 -2 

## [ 2 ,] -2 +0 

The sum of squares of the off-diagonal elements is 5. 

For 2x2 matrices we only have to compute a single plane rotation, there is no cycling. So 
let’s perform the necessary calculations. 

The matrix S is 

## C,1] [, 2] 

## [1,] +5.00 -1.00 
## [2,] -1.00 +1.25 

with smallest eigenvalue 1 and with corresponding eigenvector -0.242535625, -0.9701425001. 
Thus u = cos(2 6) is -0.242535625. 

## [1] 1 

## [1] 1 

From u we compute \^J 2(1 + u) and 2(1 — u), which are respectively 0.6154122094 and 
0.788205438. And the rotation matrix is 

## [,1] [,2] 

## [1,] +0.615412 -0.788205 
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## [2,] +0.788205 +0.615412 

The rotated matrices are 

## [, 1 ] [, 2 ] 

## [1,] +0.029857 +0.242536 
## [2,] +0.242536 +1.970143 

## [, 1 ] [, 2 ] 

## [1,] +0.757464 -0.970143 
## [2,] -0.970143 +1.242536 

## [, 1 ] [, 2 ] 

## [1,] -1.561553 +0.000000 
## [2,] +0.000000 +2.561553 

with a sum of squares of off-diagonal elements of 1 and a sum of squares of diagonal elements 
of 15. 


4 Examples 

4.1 Three Two by Two Matrices 

This is the same example we used as an illustration earlier. 

a <- c(l,- 1 , 1 , 2 , 0 , 0 , 1 ,- 2 , 0 ) 
h <- sjacobi (a, 2, 3) 


After 2 cycles (where the second one merely concludes there is convergence) we find the 
rotation 

## [, 1 ] [, 2 ] 

## [1,] 0.7882054380 0.6154122094 

## [2,] -0.6154122094 0.7882054380 

The sum of squares of the diagona elements function increased from 7 to 15. The rotated 
matrices are 


## 


[, 1 ] 

[, 2 ] 

## 

[1J 

+1.970143 

-0.242536 

## 

[ 2 ,] 

-0.242536 

+0.029857 

## 


[, 1 ] 

[, 2 ] 

## 

[1J 

+1.242536 

+0.970143 

## 

[ 2 ,] 

+0.970143 

+0.757464 

## 


[, 1 ] 

[, 2 ] 

## 

[1J 

+2.561553 

-0.000000 

## 

[ 2 ,] 

-0.000000 

-1.561553 
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The results are the same as those computed “by hand”. 


4.2 One Matrix 

The algorithm can obviously also be applied if m — 1, in which case it becomes the cyclic 
version of the Jacobi method. Minimum loss should be zero. 


## 


[, 1 ] 

[, 2 ] 

[, 3] 

[,4] 

[, 5] 

[, 6 ] 

[,7] 

[, 8 ] 

[, 9] 

[, 10 ] 

## 

[ 1 J 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

## 

[ 2 ,] 

2 

11 

12 

13 

14 

15 

16 

17 

18 

19 

## 

[3,] 

3 

12 

20 

21 

22 

23 

24 

25 

26 

27 

## 

[4,] 

4 

13 

21 

28 

29 

30 

31 

32 

33 

34 

## 

[5,] 

5 

14 

22 

29 

35 

36 

37 

38 

39 

40 

## 

[ 6 ,] 

6 

15 

23 

30 

36 

41 

42 

43 

44 

45 

## 

[7,] 

7 

16 

24 

31 

37 

42 

46 

47 

48 

49 

## 

[ 8 ,] 

8 

17 

25 

32 

38 

43 

47 

50 

51 

52 

## 

[9,] 

9 

18 

26 

33 

39 

44 

48 

51 

53 

54 

## 

[ 10 ,] 

10 

19 

27 

34 

40 

45 

49 

52 

54 

55 


h <- sjacobi (a, 10, 1) 

In 26 cycles we reduce the sum of squares of the off-diagonal elements from 84636 to 
0.0000000003. The diagonal of the rotated matrix is 


## 

[1] 

1.0699214091 

12.1639813624 

314.7797170547 

2.1774756456 

## 

[5] 

2.8050481734 

0.5991942823 

6.6137980129 

-1.8824366513 

## 

[9] 

0.1409608363 

1.5323398746 




which can be compared with the eigenvalues computed by R 


## 

[1] 

314.7797170547 

12.1639813624 

6.6137980129 

2.8050481734 

## 

[5] 

2.1774756456 

1.5323398746 

1.0699214091 

0.5991942823 

## 

[9] 

0.1409608363 

-1.8824366513 




Clearly the eigenvalues computed by sjacobi and by eigen are the same, except for order. 
And so will be the eigenvectors. 

4.3 General Case 

For our last example we construct 4 commuting matrices of order 4. Minimum loss should be 
zero, again. 

set.seed(12345) 

cl <- crossprod (matrix(rnorm(40) , 10, 4)) 
ee <- eigen(cl) $vectors 

c2 <- tcrossprod (ee °/ 0 *°/o diag(rnorm(4) ) , ee) 

c3 <- tcrossprod (ee °/ 0 *°/o diag(rnorm(4) ) , ee) 

c4 <- tcrossprod (ee diag(rnorm(4) ), ee) 
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sc <- sum(cl~2) + sum (c2~2) + sum (c3~2) + sum (c4~2) 

a <- c (cl [outer (1:4, 1: 4, ">=")] ,c2 [outer (1:4,1:4, ">=")] ,c3 [outer (1:4,1:4, ">=")] ,c4 [outer 
h <- sjacobi (a, 4, 4) 

After 4 iterations we have reduced loss from 227.4632340211 to 0.0000000000. Close enough 
to zero. 


5 Code 

5.1 R code 

5.1.1 sjacobi.R 

dyn.load(" sjacobi .so") 

sjacobi <- 
function (a, 
n, 
m, 

eps = le-15, 
itmax = 1000, 
vectors = TRUE, 
verbose = FALSE) { 

h <- .C( 

"sjacobi" , 
n = as.integer (n), 
m = as.integer (m), 
a = as.double (a), 
k = as.double (rep(0, n * n)), 
fstart = as.double (0), 
ffinal = as.double (0), 
itel = as.integer (0), 
itmax = as.integer (itmax), 
eps = as.double (eps), 
verbose = as.integer (verbose), 
vectors = as.integer (vectors) 

) 

return(list ( 
astart = a, 
afinal = h$a, 
k = h$k, 

fstart = h$f start, 
ffinal = h$f final, 
itel = h$itel 
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)) 

} 

triangle2matrix <- function (x) { 
m <- length (x) 

n <- round ((sqrt (1 + 8 * m) - 1) / 2) 
h <- 

.CC'trimat", as.integer (n), as.double (x), as.double (rep (0, n * n))) 
return (matrix (h[[3]], n, n)) 

> 

matrix2triangle <- function (x) { 
n <- dim(x) [1] 
m <- n * (n + 1) / 2 
h <- 

. CC'mattri", as.integer (n), as.double (x), as.double (rep (0, m))) 
return (h[ [3] ]) 

> 

trianglePrint <- function (x, 

width = 6, 
precision = 4) { 

m <- length (x) 

n <- round ((sqrt (1 + 8 * m) - 1) / 2) 
h <- 

. C( "pritru" , 

as.integer (n), 
as.integer (width), 
as.integer (precision), 
as.double (x)) 

} 

matrixPrint <- function (x, 

width = 6, 
precision = 4) { 

n <- nrow (x) 
m <- ncol (x) 
h <- 
• C( 

"primat" , 

as.integer (n), 

as.integer (m), 

as.integer (width), 

as.integer (precision), 
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as.double (x) 


> 


) 


5.2 C Code 

5.2.1 sjacobi.h 

#ifndef SJACOBI_H 
#define SJACOBI_H 

#include <math.h> 
#include <stdbool.h> 
#include <stdio.h> 

#include <stdlib.h> 


static 

inline 

int 

VINDEX (const 

int) 






static 

inline 

int 

MINDEX (const 

int , 

const 

int, const 

int) ; 



static 

inline 

int 

SINDEX (const 

int , 

const 

int, const 

int) ; 



static 

inline 

int 

TINDEX (const 

int , 

const 

int, const 

int) ; 



static 

inline 

int 

AINDEX (const 

int , 

const 

int, const 

int, const 

int, const 

int) ; 

static 

inline 

int 

UINDEX (const 

int , 

const 

int, const 

int, const 

int, const 

int) ; 


static 

static 

static 

static 


inline double SQUARE (const double); 
inline double THIRD(const double); 
inline double FOURTH (const double) ; 
inline double FIFTH (const double) ; 


static 

static 

static 

static 


inline double MAX(const double, const double); 
inline double MIN(const double, const double); 
inline int IMIN (const int, const int) ; 
inline int IMAX (const int, const int); 


static inline int VINDEX (const int i) { return i - 1; } 

static inline int MINDEX(const int i, const int j, const int n) { 
return (i - 1) + (j - 1) * n; 

} 


static inline int AINDEX(const int i, const int j, const int k, const int n, 

const int m) { 

return (i - 1) + (j - 1) * n + (k - 1) * n * m; 

> 


10 


static inline int SINDEX (const int i, const int j, const int n) { 

return ((j - 1) * n) - (j * (j - 1) / 2) + (i - j) - 1 ; 

> 

static inline int TINDEX(const int i, const int j, const int n) { 

return ((j - 1) * n) - ((j - 1) * (j - 2) / 2) + (i - (j - 1)) - 1; 

> 

static inline int UINDEX(const int i, const int j, const int k, const int n, 

const int m) { 

return ((k - 1) * n * (n + 1) / 2) + ((j - 1) * n) - ((j - 1) * (j - 2) / 2) + 
(i - (j - 1)) - 1; 

} 

static inline double SQUARE (const double x) { return x * x; }- 
static inline double THIRD(const double x) { return x * x * x; } 
static inline double FOURTH (const double x) { return x*x*x*x; } 
static inline double FIFTH(const double x) { return x*x*x*x*x; } 

static inline double MAX (const double x, const double y) { 
return (x > y) ? x : y; 

} 

static inline double MIN (const double x, const double y) { 
return (x < y) ? x : y; 

} 

static inline int IMAX (const int x, const int y) { return (x > y) ? x : y; } 

static inline int IMIN (const int x, const int y) { return (x < y) ? x : y; } 

void sjacobi( const int *, const int *, double *, double double *, double *, 

int const int *, const double const bool *, const bool *); 


void primat (const 

int 

*, 

const 

int *, 

const int 

*, 

const 

int *, 

const 

double 

*); 

void pritru(const 

int 

*, 

const 

int *, 

const int 

*, 

const 

double 

*); 



void prisru(const 

int 

*, 

const 

int *, 

const int 

*, 

const 

int *, 

const 

double 

*); 

void trimat (const 

int 

*, 

const 

double 

*, double 

*); 






void mattri (const 

int 

*, 

const 

double 

*, double 

*); 







#endif /* SJACOBI_H */ 
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5.2.2 sjacobi.c 


#include "sjacobi.h" 

/* 

int main(void) { 

int n = 10, m = 6, it = 0, Umax = 100, w = 10, p = 6; 
bool verbose = true, vectors = true; 

double a[330], evec[100], oldi [10], old][10], fini = 0, f = 0, eps = le-15; 
for (int i = 1; i <= 330; i++) 
a[VINDEX(i)] = (double)i; 

(void)sjacobi(&n, &m, a, evec, &fini, &f, &it, oldi, old), &itmax, &eps, 
&verbose, &vectors); 

(void)primat(&n, &n, &w, &p, evec); 

(void)prisru(&n, &m, &w, &p, a); 

} 

*/ 


int main(void) { 

int n = 2, m = 3, it = 0, Umax = 100; 
bool verbose = true, vectors = true; 

double a[9] = {1.0, -1.0, 1.0, 2.0, 0.0, 0.0, 1.0, -2.0, 0.0}; 
double evec[f] = {0.0, 0.0, 0.0, 0.0}; 

double oldi[2] = {0.0, 0.0}, old][2] = {0.0, 0.0}, fini = 0, f = 0, 
eps = le-15; 

(void)sjacobi(&n, &m, a, evec, &fini, &f, &it, oldi, oldj, &itmax, &eps, 
&verbose, &vectors); 


} 

*/ 


int main(void) { 

int n = 10, m = 1, it = 0, Umax = 100, w = 10, p = 6; 
bool verbose = true, vectors = true; 

double a [55], oldi [10], evec [100], oldj [10], fini = 0, f = 0, eps = le-16; 
for (int i = 1; i <= 55; i++) { 
a[VINDEX(i)] = (double)i; 

} 

(void)sjacobi(&n, &m, a, evec, &fini, &f, &it, oldi, oldj, &itmax, &eps, 
&verbose, &vectors); 

(void)pritru(&n, &w, &p, a); 

(void)primat(&n, &n, &w, &p, evec); 
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*/ 

void sjacobi( const int *nn, const int *mm, double *a, double *evec, 
double *fini, double *f, int *it, const int *itmax, 
const double *eps, const bool *verbose, const bool *vectors) { 
int n = *nn, m = *mm, itel = 1; 

double d = 0.0, cost = 0.0, sint = 0.0, u = 0.0, p = 0.0, q = 0.0, r = 0.0, 
piil = 0.0, pijl = 0.0, lbd = 0.0, dd = 0.0, pp = 0.0, dp = 0.0; 
double fold = 0.0, fnew = 0.0, oldi = 0.0, oldj = 0.0; 
if (^vectors) { 

for (int i = 1; i <= n; i++) { 

for (int j = 1; j <= n; j++) { 

evec[MINDEX(i, j, n)] = (i == j) ? 1.0 : 0.0; 

} 

> 

} 

for (int k = 1; k <= m; k++) { 

for (int i = 1; i <= n; i++) { 

fold += SQUARE(a[UINDEX(i, i, k, n, m)]); 

> 

} 

*fini = fold; 
while (true) { 

for (int j = 1; j <= n - 1; j++) { 
for (int i = j + 1; i <= n; i++) { 
dd = 0.0, pp = 0.0, dp = 0.0; 
for (int k = 1; k <= m; k++) { 
p = a[UINDEX(i, j, k, n, m)]; 
q = a[UINDEX(i, i, k, n, m)]; 
r = a[UINDEX(j, j, k, n, m)]; 
d = (q - r) / 2.0; 
dd += SQUARE(d); 
pp += SQUARE(p); 
dp += p * d; 

} 

lbd = ((dd + pp) - sqrt(SQUARE(dd - pp) + 4.0 * SQUARE(dp))) / 2.0; 
u = dp / sqrt(SQUARE(dp) + SQUARE(lbd - pp)); 
if ((fabs(dp) < le-15) && (pp <= dd)) { 
continue ; 

} 

cost = sqrt((l + u) / 2); 
sint = -sqrt((l - u) / 2); 
if (^vectors) { 

for (int 1=1; 1 <= n; 1++) { 
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piil = evec[MINDEX(1, i, n)]; 
pijl = evec [MINDEX(1, j, n)]; 

evec[MINDEX(1, i, n)] = cost * piil - sint * pijl; 
evec [MINDEX(1, j, n)] = sint * piil + cost * pijl; 

> 

} 

for (int k = 1; k <= m; k++) { 
p = a [UINDEXd, j, k, n, m)]; 

q = a [UINDEXd, i, k, n, m) ] ; 

r = a[UINDEX(j, j, k, n, m)]; 

for (int 1=1; 1 <= n; 1++) { 
if ((1 == i) || (1 == j)) 
continue ; 

int il = IMINd, 1); 

int li = IMAXCi, 1); 

int jl = IMINCj, 1); 

int lj = IMAXCj, 1); 

oldi = a[UINDEX(li, il, k, n, m)]; 

oldj = a[UINDEX(1j, jl, k, n, m)]; 

a[UINDEX(li, il, k, n, m)] = cost * oldi - sint * oldj; 

a[UINDEX(lj, jl, k, n, m)] = sint * oldi + cost * oldj; 

} 

a [UINDEXd, i, k, n, m)] = 

SQUARE(cost) * q + SQUARE(sint) * r - 2 * cost * sint * p; 
a[UINDEX(j, j, k, n, m)] = 

SQUARE(sint) * q + SQUARE(cost) * r + 2 * cost * sint * p; 
a[UINDEX(i, j, k, n, m)] = 

cost * sint * (q - r) + (SQUARE(cost) - SQUARE(sint)) * p; 

} 

} 

> 

fnew = 0.0; 

for (int k = 1; k <= m; k++) { 
for (int i = 1; i <= n; i++) { 

fnew += SQUARE(a[UINDEX(i, i, k, n, m)]); 

> 

> 

if (^verbose == true) { 

printf("itel °/„4d fold %15.10f fnew °/„15. lOf \n", itel, fold, fnew); 

> 

if (((fnew - fold) < *eps) I I (itel == *itmax)) 
break; 

fold = fnew; 
itel++; 
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} 

*f = fnew; 

*it = itel; 
return; 

> 

void primat( const int *n, const int *m, const int *w, const int *p, 
const double *x) { 
for (int i = 1; i <= *n; i++) { 
for (int j = 1; j <= *m; j++) { 

printf(" %*.*f ", *w, *p, x[MINDEX(i, j, *n)]); 

> 

printf("\n"); 

} 

printf("\n\n"); 
return; 

} 

void pritru(const int *n, const int *w, const int *p, const double *x) { 
for (int i = 1; i <= *n; i++) { 
for (int j = 1; j <= i; j++) { 

printf (" %*.*f ", *w, *p, x[TINDEX(i, j, *n)]); 

> 

printf("\n"); 

} 

printf("\n\n"); 
return; 


void prisru(const int *n, const int *m, const int *w, const int *p, 
const double *x) { 
for (int k = 1; k <= *m; k++) { 
for (int i = 1; i <= *n; i++) { 
for (int j = 1; j <= i; j++) { 

printf (" °/ 0 *.*f ", *w, *p, x[UINDEX(i, j, k, *n, *m)]); 

} 

printf("\n"); 

> 

printf("\n\n"); 

} 

return; 


void trimat (const int *n, const double *x, double *y) { 
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int nn = *n; 

for (int i = 1; i <= nn; i++) { 

for (int j = 1; j <= nn; j++) { 

y[MINDEX(i, j, nn)] = 

(i >= j) ? x[TINDEX(i, j, nn)] : x[TINDEX(j, i, nn)]; 

> 

} 

return; 

} 

void mattri( const int *n, const double *x, double *y) { 
int nn = *n; 

for (int j = 1; j <= nn; j++) { 

for (int i = j; i <= nn; i++) { 

y[TINDEX(i, j, nn)] = x[MINDEX(i, j, nn)]; 

> 

} 

return; 

> 
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